/*******************************************************************************
*
*   Copyright ©2009 Alexandru Marin
*
*   This file is part of nwave/pMRI C library.
*
*   This project started as an MSc degree project during a
*   stage at the Paris-Est university, under joint supervision of prof.
*   Bogdan Ionescu (Image Analysis and Processing Laboratory, "Politehica"
*   University of Bucharest), Caroline Chaux and Jean-Christophe Pesquet
*   (Laboratoire d'Informatique, Université Paris-Est, Institut Gaspard Monge).
*
*   nwave/pMRI C libray is free software: you can redistribute it and/or modify
*   it under the terms of the GNU General Public License as published by
*   the Free Software Foundation, either version 2 of the License, or
*   (at your option) any later version.
*
*   nwave/pMRI C libray is distributed in the hope that it will be useful,
*   but WITHOUT ANY WARRANTY; without even the implied warranty of
*   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*   GNU General Public License for more details.
*
*   You should have received a copy of the GNU General Public License
*   along with nwave/pMRI C library.  If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/

#include "mri_types.h"
#include "mri_utils.h"
#include <assert.h>
#include <stdlib.h>

mri_image_t* mri_image_read(FILE* f) {

	mri_image_t* result = (mri_image_t*)malloc(sizeof(mri_image_t));
	assert(result != NULL);

	mri_matrix_descriptor_t* desc = (mri_matrix_descriptor_t*)malloc(sizeof(mri_matrix_descriptor_t));
	assert(desc != NULL);

	result->desc = desc;

	result->data = mri_read_complex_matrix(f, desc);
	assert(desc->ndims > 1 && desc->ndims < 5);
	result->width = desc->dims[0];
	result->height = desc->dims[1];
	if(desc->ndims > 2)
		result->nslices = desc->dims[2];
	else
		result->nslices = 1;
	if(desc->ndims > 3)
		result->nsamples = desc->dims[3];
	else
		result->nsamples = 1;

	return result;
}

mri_image_t* mri_image_load(char* name) {

	char* filename = mri_strcat(3, mri_get_path(), name, ".csv");

	FILE* f = fopen(filename, "r");
	assert(f != NULL);

	mri_image_t* result = mri_image_read(f);

	free(filename);
	fclose(f);

	return result;
}

void mri_image_write(mri_image_t* image, FILE* f) {

	assert(image!=NULL);

	int dims[] = {image->width, image->height, image->nslices, image->nsamples};
	mri_matrix_descriptor_t desc;
	desc.ndims = sizeof(dims)/sizeof(dims[0]);
	desc.dims = dims;

	mri_write_complex_matrix(f, image->data, &desc);
}

void mri_image_save(mri_image_t* image, char* name) {

	char* filename = (char*)mri_strcat(3, mri_get_path(), name, ".csv");

	FILE* f = fopen(filename, "w");
	assert(f != NULL);

	mri_image_write(image, f);

	free(filename);
	fclose(f);
}

void mri_images_write(int nimages, mri_image_t** images, FILE* f) {

	int dims[] = {images[0]->width, images[0]->height, images[0]->nslices, images[0]->nsamples, nimages};
	mri_matrix_descriptor_t desc;
	desc.ndims = sizeof(dims)/sizeof(dims[0]);
	desc.dims = dims;

	int i, j, totalsize=1;

	// Compute total size of one mri image
	for (i = 0; i < desc.ndims-1; i++)
		totalsize *= desc.dims[i];

	// Write header. First number of dimensions
	fprintf(f, "%d", desc.ndims);
	fputc(',', f);
	// Then all dimensions
	for (i = 0; i < desc.ndims; i++) {
		fprintf(f, "%d", desc.dims[i]);
		fputc(',', f);
	}

	// Write actual data
	for(j=0;j<nimages;j++)
		for (i = 0; i < totalsize; i++) {
			mri_write_complex(f, images[j]->data[i]);
			if (j != nimages - 1 || i != totalsize - 1)
				fputc(',', f);
		}

}

void mri_images_save(int nimages, mri_image_t** images, char* name) {

	char* filename = mri_strcat(3, mri_get_path(), name, ".csv");

	FILE* f = fopen(filename, "w");
	assert(f != NULL);

	mri_images_write(nimages, images, f);

	free(filename);
	fclose(f);
}

mri_image_t* mri_image_alloc(int height, int width, int nslices, int nsamples) {

	mri_image_t* result = (mri_image_t*) malloc(sizeof(mri_image_t));
	assert(result != NULL);

	result->height = height;
	result->width = width;
	result->nslices = nslices;
	result->nsamples = nsamples;

	result->data = (complex float*) malloc(MRI_IMAGE_SIZE(result) * sizeof(complex float));
	assert(result->data != NULL);

	return result;
}

mri_image_t* mri_image_gen_synth(int height, int width, int nslices, int nsamples, complex float fillvalue) {

	int i, j, k, l;

	mri_image_t* result = mri_image_alloc(height, width, nslices, nsamples);

	// Fill with data
	for (i = 0; i < result->height; i++)
		for (j = 0; j < result->width; j++)
			for (k = 0; k < result->nslices; k++)
				for (l = 0; l < result->nsamples; l++)
					MRI_IMAGE(result, i, j, k, l) = fillvalue;

	return result;
}

void mri_image_free(mri_image_t* image) {
	free(image->data);
	free(image);
}


mri_sens_t* mri_sens_alloc(int height, int width, int ncoils, int nslices) {

	mri_sens_t* result = (mri_sens_t*)malloc(sizeof(mri_sens_t));
	assert(result != NULL);

	result->height = height;
	result->width = width;
	result->ncoils = ncoils;
	result->nslices = nslices;

	result->data = (complex float*)malloc(MRI_SENS_SIZE(result)*sizeof(complex float));
	assert(result->data != NULL);

	return result;
}

mri_sens_t* mri_sens_gen_synth(int height, int width, int ncoils, int nslices, complex float fillvalue) {

	int i,j,k,l;

	mri_sens_t* result = mri_sens_alloc(height, width, ncoils, nslices);

	for(i=0;i<height;i++)
		for(j=0;j<width;j++)
			for(k=0;k<ncoils;k++)
				for(l=0;l<nslices;l++)
					MRI_SENS(result, i, j, k, l) = fillvalue;

	return result;
}

mri_sens_t* mri_sens_read(FILE* f) {

	mri_matrix_descriptor_t desc;

	mri_sens_t* result = (mri_sens_t*)malloc(sizeof(mri_sens_t));
	assert(result != NULL);

	result->data = mri_read_complex_matrix(f, &desc);
	assert(desc.ndims > 2 && desc.ndims < 5);
	result->width = desc.dims[0];
	result->height = desc.dims[1];
	if(desc.ndims > 3) {
		result->nslices = desc.dims[2];
		result->ncoils = desc.dims[3];
	} else {
		result->nslices = 1;
		result->ncoils = desc.dims[2];
	}
	return result;
}

mri_sens_t* mri_sens_load(char* name) {

	char* filename = mri_strcat(3, mri_get_path(), name, ".csv");

	FILE* f = fopen(filename, "r");
	assert(f != NULL);

	mri_sens_t* result = mri_sens_read(f);

	free(filename);
	fclose(f);

	return result;
}

void mri_sens_free(mri_sens_t* sens) {
	free(sens->data);
	free(sens);
}
